/*****************************************************************************
 * $LastChangedDate: 2010-10-23 18:22:39 -0400 (Sat, 23 Oct 2010) $
 * @file
 * @author  Jim E. Brooks  http://www.palomino3d.org
 * @brief   Math functions for computer graphics.
 * @remarks - A few functions are named specially to avoid ambiguity with standard functions.
 *          - To speed or improve FP execution:
 *            - use multiplication instead of division
 *            - power of 2 operations involve no loss of accuracy
 *            - math::PI is preferred over M_PI etc from math.h.
 *              math.h has confusing #defines such as M_2_PI
 *              which means 2/PI (not 2*PI).
 *          - NaN are caused by division and some standard C math functions.
 *          - All modern FPUs are based on IEEE 754.
 *            A single-precision 32-bit FP is limited to 24 bits (8 decimal digits) of precision.
 *          - C++ <cmath> is better than C <math.h> because C++ functions are overloaded.
 * REMINDERS: - Consider using MATH_EXCEPTION() instead of ASSERT().
 *//*
 * LEGAL:   COPYRIGHT (C) 2004 JIM E. BROOKS
 *          THIS SOURCE CODE IS RELEASED UNDER THE TERMS
 *          OF THE GNU GENERAL PUBLIC LICENSE VERSION 2 (GPL 2).
 *****************************************************************************/

#ifndef MATH_FUNCS_HH
#define MATH_FUNCS_HH 1

#if _MSC_VER
#include <float.h>  // _finite()
#endif // _MSC_VER

#include "math/funcs_debug.hh"

namespace math {

////////////////////////////////////////////////////////////////////////////////
//////////////////////////  D E F I N I T I O N S  /////////////////////////////
////////////////////////////////////////////////////////////////////////////////

namespace defs {

const fpx FP_DELTA_TOLERANCE = 0.001;   ///< used for comparing fp values

}  // defs

////////////////////////////////////////////////////////////////////////////////
///////////////////////////  C O M P A R I S O N  //////////////////////////////
////////////////////////////////////////////////////////////////////////////////

/*****************************************************************************
 * IfNan() returns true if a Not-A-Number or Infinity.
 * C99 isfinite() is true if neither NaN nor infinity.
 * Technically, divide-by-zero results in Infinity, not NaN,
 * but IfNan() is used to check divide errors.
 * A NaN never equals any number -- not even itself.
 * !(n==n) is unreliable because smart compilers will optimize as always false (gcc will).
 *****************************************************************************/
#if _MSC_VER
#   define IfNan(N) (UX(!_finite(N)))
#else
#   define IfNan(N) (UX(!std::isfinite(N)))
#endif

/*****************************************************************************
 * @return True if two fp values are approximately equal.
 * NOTE: Will malfunction if arg is NaN/Infinity.
 *****************************************************************************/
template<typename FP> bool FP_EQ( FP a, FP b )
{
ASSERT( not (IfNan(a) or IfNan(b)) );

    return std::abs(a-b) < math::defs::FP_DELTA_TOLERANCE;
}

/*****************************************************************************
 * @return True if one fp value is less than or approximately equal to another.
 * NOTE: Will malfunction if arg is NaN/Infinity.
 *****************************************************************************/
template<typename FP> bool FP_LE( FP a, FP b )
{
ASSERT( not (IfNan(a) or IfNan(b)) );

    return (a - math::defs::FP_DELTA_TOLERANCE <= b)
        or (a + math::defs::FP_DELTA_TOLERANCE <= b);
}

/*****************************************************************************
 * @return True if one fp value is greater than or approximately equal to another.
 * NOTE: Will malfunction if arg is NaN/Infinity.
 *****************************************************************************/
template<typename FP> bool FP_GE( FP a, FP b )
{
ASSERT( not (IfNan(a) or IfNan(b)) );

    return (a - math::defs::FP_DELTA_TOLERANCE >= b)
        or (a + math::defs::FP_DELTA_TOLERANCE >= b);
}

/*****************************************************************************
 * @return True If two numbers have the same sign.
 *****************************************************************************/
template<typename T1,typename T2>
bool
IfSameSign( const T1 a, const T2 b )
{
  //return (a * b) >= 0;  // wrong because if either is zero: -n * 0 = 0

    if ( (a >= 0) and (b >= 0) )  // both positive?
        return true;
    else if ( (a < 0) and (b < 0) )  // both negative
        return true;
    else
        return false;  // different signs
}

////////////////////////////////////////////////////////////////////////////////
///////////////////////////////  B A S I C   ///////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

/*****************************************************************************
 * Square root function.
 *****************************************************************************/
template<typename FP>
FP
SquareRoot( FP f )
{
MATH_EXCEPTION( f >= 0.0f and not IfNan(f), "math error: sqrt(-n|NaN) undefined" );
    return std::sqrt( f );  // C++ <cmath>
}

////////////////////////////////////////////////////////////////////////////////
////////////////////////////  M I S C   M A T H   //////////////////////////////
////////////////////////////////////////////////////////////////////////////////

/*****************************************************************************
 * x^2
 *****************************************************************************/
template<typename T>
T
Square( T x )
{
    return x * x;
}

/*****************************************************************************
 * Average.
 *****************************************************************************/
template<typename T>
T
Average( T x, T y )
{
    return (x + y) / 2;
}

/*****************************************************************************
 * Force a number to be within a range.
 *****************************************************************************/
template<typename T>
T
Range( T val, T lo, T hi )
{
    if ( val < lo )
        return lo;
    else if ( val > hi )
        return hi;
    else
        return val;
}

/*****************************************************************************
 * @return True if two lines overlap.
 *****************************************************************************/
template<typename T>
T
IfInRange( T a0, T a1, T b0, T b1 )
{
    // Prove the reverse: test if out-of-range.
         if ( a0 > b1 )
        return false;
    else if ( a1 < b0 )
        return false;
    else if ( b0 > a1 )
        return false;
    else if ( b1 < a0 )
        return false;
    else
        return true;  // within range
}

/*****************************************************************************
 * Compute a fraction for a number in a range.
 * Number is allowed to be outside the range (result will be clamped to 0 or 1).
 * Examples:
 * FractionInRange(100,100,200) = 0.0
 * FractionInRange(150,100,200) = 0.5
 * FractionInRange(200,100,200) = 1.0
 *****************************************************************************/
INLINE fp
FractionInRange( fp n, fp low, fp high )
{
// Negative range undefined.
MATH_EXCEPTION( (low < high) and (low >= 0 and high > 0),
                "math error: FractionInRange()" );

    if ( n <= low )
    {
        return 0.0f;
    }
    else if ( n >= high )
    {
        return 1.0f;
    }
    else  // inside range
    {
        return (n - low) / (high - low);
    }
}

/*****************************************************************************
 * Map a number from one range into another.
 *****************************************************************************/
INLINE fp
Remap( fp value, fp rangeOld, fp rangeNew )
{
    const fp n = value / rangeOld * rangeNew;
MATH_EXCEPTION( not IfNan(n), "math error: Remap()" );
    return n;
}

/*****************************************************************************
 * Truncate a float at an interval.
 * Eg, Truncate(24,10) returns 20.
 *****************************************************************************/
INLINE fp
Truncate( fp val, fp interval )
{
    float integral;  // float, not fp, as address is passed to modff()
    std::modf( val/interval, &integral );
    return integral * interval;
}

/*****************************************************************************
 * Modulo that handles negative values by subtracting from limit.
 * For example, ModuloAbs(deg,360.0f) can be used to accomodate a function
 * that requires degrees to be in range {0,..,360}
 * ModuloAbs( -1.0,  360.0 ) = 359.0
 * ModuloAbs( -90.0, 360.0 ) = 270.0
 *****************************************************************************/
template<typename T>
T
ModuloAbs( const T n, const T limit )
{
    const T m = std::fmod( n, limit );
    if ( m >= 0 )
        return m;
    else  // m < 0
        return limit + m;  // effectively subtraction
}

} // namespace math

#endif // MATH_FUNCS_HH
